NPS-PH-92-008 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


A  Self-affine  Multi-fractal  Wave/ 

Turbulence  Discrimination  Method  Using  Data 

from  Single  Point  Fast  Response  Sensors 

in  a  Nocturnal  Atmospheric  Boundary  Layer 


R.F 


by 
Kamada  and  A.J.  DeCaria 
April  1992 


Report  for  Period 
October  1990  -  September  1991 


FEDDOCS 

D  208.14/2 
NPS-PH-92-008 


Approved  for  public  release;  distribution  unlimited 

Prepared  for:  U.S.  Air  Force,  Space  Division, 

Los  Angeles  AFB,  California  90009 


fa 

"^RCT.CWJFO"0**"**^     NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 

Rear  Admiral    R.  W.  West  "•   shu11 

Superintendent  Provost 

The  work  reported  herein  was  prepared  for  and  funded  by  U.S.  Air 
Force  Space  Systems  Division,  Los  Angeles,  CA,  90009. 

Reproduction  of  all  or  part  of  this  document  is  authorized. 
This  report  was  prepared  by 


UNCLASSIFIED 


ECURITY   CLASSiF'CATIQN   Oc    TH'S    PAGE 


***m-  ^'SiaRADUAtWtCHOOI. 
MOWTERFY  CAHFORNrAa3343W»» 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No  0704-0188 


la    REPORT  SECURITY  CLASS  FiCATION 

UNCLASSIFIED 


lb    RESTRICTIVE  MAR^NCS 


l.'a    SECURiTY  CLASSIFICATION  AUTHORITY 


?b    DECLASSIFICATION/ DOWNGRADING  SCHEDULE 


3    Distribution/ ava.labi.my  of  report 
Approved   for  public   release; 
distribution    is   unlimited 


J    PERFORMING  ORGANIZATION  REPORT  NUMBERS) 

NPS-PH-92-008 


5    MONITORING  ORGANIZATION  REPORT  NuM3ER;S) 


5a    NAME  OF  PERFORMING  ORGANIZATION 

Naval   Postgraduate   School 


6b    OFFiCE   SYMBOL 
(If  applicable) 


la    NAME  OF  MONITORING  ORGANiZA"  ON 

USAF  Space  Systems  Division 


6c.  ADDRESS  (City,  State,  and  ZIP  Code) 


Monterey,    CA     93943-5000 


7b    ADDRESS  (Cry,  Stare,  and  ZIP  Code) 


Los   Angeles  AFB ,    CA   90009 


3a    NAME  OF  FUNDING  .'SPONSORING 
ORGANIZATION 

USAF  Space  Systems  Division 


8d    OFFICE  SYMBOL 
(If  applicable) 

USAF/SSD 


9    PROCUREMENT  INSTRUMENT   IDENTIFICATION   NUMBER 


MPIR  FY76169100412 


8c    ADDRESS  (City,  State,  and  ZIP  Code) 


Los  Angeles  AFB,    CA      90009 


10  source  of  funding  numbers 


PROGRAM 
ELEMENT  no 


2R0>ECT 

no 


task 

NO 


*VOF<    UNiT 
ACCESSION  NO 


11  TITLE  (Include  Security  Classification) 

A  Self-affine  Multi-fractal  Wave/Turbulence  Discrimination  Method  Using  Data  from  Single 
Point  Fast  Response  Sensors  in  a  Nocturnal  Atmospheric  Boundary  Layer 

12  PERSONAL  AUThOR(S) 

R.F.    Kamada 


13a    TYPE  OF  REPORT 

Technical  Report 


3b    TIME  COVERED 

from  10-90      to  9-91 


4    DATE  OF  REPORT    (Year,  Month  Day)     |15    PAGE  COUNT 

April    10,    1992 j        45 


16    SUPPLEMENTARY   NOTATION 


17 


COS* 


CODES 


FIELD 


GROUF 


SUB-GROUP 


18    SUBjECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Stable  boundary  layer,  self-affine  fractals,  turbulence 
intermittency ,  coherent  structures,  wave  trains,  wave- 
breaking 


19    ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

We  present  DA,  a  self-affine,  multi-fractal  which  may  become  the 
first  routine  wave/turbulence  discriminant  for  time  series  data. 
Using  nocturnal  atmospheric  data,  we  show  the  advantages  of  DA  over 
self-similar  fractals  and  standard  turbulence  measures  such  as 
FFTs,  Richardson  number,  Brunt-Vaisala  frequency,  buoyancy  length 
scale,  variances,  turbulent  kinetic  energy,  and  phase  averaging. 
DA  also  shows  promise  in  resolving  "wave-break"  events.  Since  it 
uses  local  basis  functions,  DA  may  be  an  ideal  tool  to  detect 
intermittent  turbulence,  coherent  structures,  and  discrete  wave 
trains  in  general.   DA  may  also  be  a  measure  of  chaos  in  general. 


20    DISTRIBUTION  'AVAILABILITY  OF  ABSTRACT 
^UNCLASSIFIED  UNLIMTED       D   SAME   AS  RPT  D  DTIC  USERS 


21     ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


12a    NAME  OF  RESPONSIBLE   INDIVIDUAL 

R.F.    Kamada 


22b  TELEPHONE  (Include  Area  Code) 

(408)    646-2674 


22c    OFFICE  SYMBOL 
PH 


DDForm  1473,   JUN  86 


Previous  editions  are  obsolete 
S/N    0102-LF-014-6603 


SECURITY  CLASSIFICATION  OF  THiS  PAGE 

UNCLASSIFIED 


I.  INTRODUCTION   1 

II.  THEORY 3 

A.  Current  ways  to  distinguish  waves  from  turbulence  3 

B.  Basic  concepts  of  self-similar  fractal  dimension.  .  3 

C.  Time  series  applications  9 

D.  Relevant  concepts  of  gravity  waves  and  turbulence.  .  10 

III.  METHODS 12 

A.  Data 12 

B.  Analysis  Measures  14 

IV.  RESULTS  AND  DISCUSSION 20 

A.  Comparisons  with  local  measures 20 

B.  Comparisons  with  bulk  measures 33 

IV.   SUMMARY  AND  CONCLUSIONS 39 

REFERENCES 4  1 

DISTRIBUTION  LIST 43 


I.  INTRODUCTION 

The  second,  third,  and  fourth  moments  of  fluctuation  velocity, 
temperature,  and  passive  tracer  concentrations,  and  their  cross 
correlations  have  generally  been  used  to  estimate  the  magnitude  of 
dispersive  atmospheric  transport.  This  applies  to  the  transport  of 
momentum,  energy,  humidity,  and  contaminants  such  as  pollutants. 
To  avoid  overestimates  of  transport,  we  must  also  be  able  to 
distinguish  between  wave  induced  and  turbulent  contributions  to 
these  fluctuation  quantities,  since  periodic  wave  motion  is  much 
less  dispersive  than  aperiodic  turbulence.  This  applies 
particularly  in  transitionally  turbulent  regimes  and  regions  such 
as  the  nocturnal  atmospheric  boundary  layer  and  convective 
entrainment  zone,  where  the  mean  flow  is  punctuated  by  both 
intermittent  turbulence  and  waves.  However,  current  algorithms  for 
separating  waves  from  turbulence  are  either  unreliable,  not 
operationally  useful,  or  both.  a  facile  wave/turbulence 
discriminant  has  not  been  demonstrated  in  the  literature. 

Here  we  study  DA,  a  self-affine,  multi-fractal  which  may  be  the 
first  operationally  useful  wave/turbulence  discriminant  for  typical 
time  series  data.  We  discuss  and  demonstrate  its  advantages  over 
a  self-similar  fractal,  Fourier  transforms,  and  other  standard 
turbulence  measures,  such  as  Richardson  number  (R,)  ,  Brunt-Vaisala 
frequency  (BVF)  ,  buoyancy  length  scale  (1B)  ,  variances,  and 
turbulent  kinetic  energy  (TKE) .  We  also  discuss  this  new  method's 
operational  advantages  over  the  phase  averaging  technique.  For  the 
comparisons  we  use  nocturnal  data  from  the  300  meter  tower  at  the 
Boulder  Atmospheric  Observatory.  We  also  suggest  that  DA  be  tested 
as  a  general  measure  of  the  degree  of  chaos  in  systems  ranging  from 
waves  and  limit  cycles  to  standard  Lorenz,  Henon,  and  Poincare 
maps,  to  transitional  turbulence,  and  turbulence  displaying  an 
extended  inertial  subrange. 

Fractal  self-similarity  exists  when  a  geometric  feature  re-appears 
at  successively  smaller  scales,  where  each  successive  scale  is 
separated  by  a  constant  scaling  factor  (Mandelbrot,  1977).  In 
fluids  where  turbulence  lacks  rigid  geometric  structure,  we  must 
instead  test  the  statistical  properties  describing  the  ensemble 
mean  geometry  of  the  flow  structure  for  possible  self-similarity. 

Previous  studies  have  suggested  strong  potential  for  a  fractal 
self-similarity  approach  to  wave/turbulence  discrimination. 
Feigenbaum  (1980)  discovered  that  low  order  turbulence,  i.e., 
chaos,  in  viscousless  systems  possesses  an  underlying  self- 
similarity  and  hence  fractally  intermittent  character.  He  showed 
that  this  self-similarity  in  the  density  of  Lagrangian  particle 
streamlines  appears  only  when  the  attractor,  the  point  or  area 
about  which  the  particle  motion  occurs  in  real  or  phase  space 
becomes  unstable,  i.e.,  strange.  For  waves  or  even  limit  cycle 
motion  about  fixed  attractors,  self-similarity  does  not  occur. 
Pertinent  to  time  series  analysis,  this  self-similarity  appears  in 
real  as  well  as  phase  space. 


With  regard  to  real  viscous  fluids,  Kamada  (1988a,b)  showed  a 
length  scale  self-similarity  in  atmospheric  transitional  turbulence 
for  the  capping  inversion  above  a  convective  mixed  layer. 
Sreenivasan  and  Meneveau  (1986)  ,  and  Presad  and  Sreenivasan  (1990) 
found  that  the  interfacial  convolutions  between  turbulent  and  non- 
turbulent  regions  of  a  shear  flow  are  self-similar  and  thus  have  a 
fractal  dimension.  However,  pure  self-similarity  does  not  apply 
across  all  atmospheric  turbulence  scales  because  eddies  are 
increasingly  compressed  vertically  with  scale,  due  to  stability. 
However,  Lovejoy  and  Schertzer  (1984)  included  this  compression  by 
using  a  multi-fractal,  self-affine  stretch  modification  to  the 
standard  self-similar  fractal.  Packard  et  al.  (1980)  and  Pawelzik 
and  Schuster  (1987)  also  studied  the  fractal  dimensions  of 
attractors  inferred  from  a  time  series. 

However,  the  aim  here  is  to  find  and  utilize  the  fractal  dimension 
of  the  time  series  itself,  not  the  dimension  of  the  phase  space 
attractor  or  the  dimension  of  the  turbulence  field  in  real  space. 
Since  a  time  series  is  a  digital  sampling,  a  time  series  trace  of 
a  fractal  process  should  also  show  fractal  characteristics.  For 
example,  Carter  et  al.  (1986)  found  that  the  known  fractal 
characteristics  of  cloud  geometry  in  real  space  were  also  evident 
in  the  time  series  of  their  infrared  intensity  versus  azimuth 
viewing  angle. 

In  real  applications  fractal  analysis  has  potential  advantages  over 
standard  spectral  methods.  1)  Less  data  manipulation  is  required. 
Standard  Fourier  transform  techniques  require  that  the  data  be 
periodic  within  the  data  window  to  avoid  the  introduction  of 
fictitious  high  frequency  noise.  This  requires  tapering  the  data 
within  the  window  so  that  its  endpoints  have  the  same  value.  This 
is  not  required  for  fractal  analysis.  2)  Fractal  analysis' 
biggest  advantage  is  that  sudden  amplitude  changes,  involving  such 
phenomena  as  truncated  gravity  wave  trains,  breaking  Kelvin- 
Helmholtz  instabilities,  or  turbulent  ramp  structures  pose  no 
problems.  In  Fourier  analyses  a  linear  combination  of  wave  numbers 
much  higher  than  the  time  resolution  is  required  to  accurately 
account  for  such  near-discontinuities.  This  is  because  the  basis 
functions  for  standard  spectral  analysis:  sinusoids,  Legendre  or 
Leguerre  polynomials,  Bessel  functions,  etc.,  are  global,  i.e., 
infinite  in  length,  rather  than  discrete;  whereas  fractal  analysis 
naturally  assumes  a  discrete  Chapeau  (top  hat)  basis  function  over 
the  time  resolution  in  question.  And  this  basis  function 
inherently  spans  the  ideal  range:  from  the  chosen  span  of  the  time 
series  to  the  available  limit  of  resolution.  Therefore,  functions 
involving  frequencies  shorter  than  the  available  resolution  are  not 
required  to  portray  large  jumps  or  discontinuities.  This  suggests 
that  fractal  analysis  may  be  quite  useful  in  analyzing  discrete, 
coherent  wave  trains  and  intermittent  turbulence  structures  in 
general. 


II.  THEORY 

A.  Current  ways  to  distinguish  waves  from  turbulence 

Since  linear  waves  transport  no  heat,  the  cross-spectral  method 
assumes  for  waves  that  the  vertical  velocity  and  temperature 
exhibit  a  small  cospectrum  and  large  quadrature  spectrum,  i.e., 
temperature  lags  vertical  velocity  by  90°,  while  turbulence  shows 
a  large  cospectrum  and  small  quadrature  spectrum.  But  Finnigan 
(1988)  noted  that  near-surface  gravity  waves  are  non-linear,  so  the 
co-spectrum  is  usually  comparable  to  the  quadrature  spectrum. 

The  spectral  gap  method  assumes  a  wave/turbulence  gap  in  the 
boundary  layer  power  spectra  (Nai-Ping,  1983).  Caughey  (1977)  ties 
the  position  of  this  gap  to  the  Brunt-Vaisala  frequency  (BVF) , 
which  should  be  the  highest  frequency  of  gravity  wave  that  the 
atmosphere  will  support.  This  frequency  is  often,  but  not  always, 
lower  than  turbulence  injection  scales.  Wave  phase  velocity  also 
adds  to  the  apparent  frequency.  Thus,  spectral  gaps  in  wave  number 
may  disappear  in  the  frequency  domain  for  time  series  from  a  fixed 
sensor   (Caughey  and  Readings  ,  1975) . 

Phase  averaging  subtracts  the  wave  and  mean  components  from  the 
signal,  leaving  turbulence  as  the  remainder  (Finnigan,  ibid).  But 
one  must  assume  a  wave  exists  and  its  frequency  must  accurately 
estimated  prior  to  the  analysis.  A  surface  microbarograph  array  is 
usually  needed  to  determine  wave  frequency  reliably,  but  most  tower 
sites  lack  microbarographs .  Problems  also  arise,  if  wave 
amplitudes  change  greatly  with  time,  if  coherence  is  lost  after  a 
few  periods,  or  if  non-linear  dispersion  changes  the  wave 
frequency.  We  suggest  that  phase  averaging  is  suitable  for  some 
aspects  of  wave/turbulence  analysis,  but  only  after  wave  presence 
has  been  established,  perhaps  most  easily  by  multi-fractal 
analysis. 

The  following  discusses  how  multi-fractal  analysis  can  establish 
such  wave  presence  from  single  point  measurements,  without 
microbarographs,  or  assumptions  about  quadrature  spectra  or 
spectral  gaps. 

B.  Basic  concepts  of  self-similar  fractal  dimnesion 

For  clarity  we  outline  basic  fractal  concepts  before  describing  the 
self-affine  multi-fractal  operator  used  for  this  analysis. 
The  minimum  number  of  boxes  of  side  e    needed  to  cover  a  set  of 
points  on  a  two-dimensional  plot  scales  as 

N(c)  =   —n    ,  (II. 1) 

where  D  is  defined  as  the  capacity  dimension  of  the  set  and  F,  the 
lacunity  (Mandelbrot,  1977)  .  If  the  set  forms  a  straight  line, 
then  D  =  1  because,  if  the  box  length  is  halved,  twice  as  many 


boxes  are  needed  to  cover  the  line.  If  the  set  consists  of 
uniformly  distributed  points  in  a  plane,  then  D  =  2,  since  four 
times  as  many  boxes  are  needed. 

If  the  points  are  not  uniformly  distributed,  then  D  can  be  non- 
integer.  The  set  of  points  is  then  fractal  with  dimension  D, 
defined  in  the  limit  as  e  approaches  zero.  This  definition  is 
called  the  box  dimension,  DB,  (Mandelbrot,  1985) ,  and  is  given  as 


=  lim 

e-0 


logN(t) 

log^ 
e 


(II. 2) 


Figure  1  -  Construction  of  the  Cantor  set 


The  Cantor  set  in  fig.  1  has  a  non-integer  DB.  It  is  formed  by 
removing  the  middle  third  of  a  line  segment,  removing  the  middle 
third  of  each  remaining  segment,  and  re-iterating  to  infinity.  The 
set  is  self-similar  because  the  geometric  structure  seen  at  the 
largest  scale  re-appears  for  all  smaller  scales  removed  by 
successive  factors  of  three.  For  an  initial  unit  length,  N(l)=l 
box  of  side  e  =  1  is  needed  to  cover  the  set.  For  £=1/3,  N(l/3)=2, 
and  for  £=1/9,  N(l/9)=4.  Thus,  e  is  generally  3",  and  N(£)  is  2n, 
where  n  =  0, 1 ,  2  ,  3  ,  .  .  .  ,oo  (Grebogi  et  al.,  1987).  Noting  that  n-^» 
implies  £—0,  from  eqn.  II.  2,  DB  for  the  Cantor  set  becomes 
intermediate  between  a  point  (DB  =  0)  and  a  line  (DB  =  1)  .   That  is, 

n  =  UmlSRlH    =  l£9_2  =  0.63092...  (II.  3) 

n—  log3n    log  3 

Unlike  the  analytic  Cantor  set,  DB  must  be  evaluated  numerically 
for  most  fractal  sets.  This  is  done  by  removing  the  limit  in  II. 2 
and  reordering  terms  to  get 

log  N{c)    =  log  F  -  DBlogz    . 


(II. 4) 


Then  -DB  is  the  slope  of  the  plot  of  log  N(£)  vs.  log  £,  and  log  F 
is  the  y-intercept. 

In  practice,  real  time  series  will  contain  some  noise,  so  L(£)  will 
not  be  exact.  Thus,  multiple  values  of  £  should  be  used,  and  some 
form  of  regression  employed  to  find  the  slope  (see  section  III). 
Moreover,  the  data  may  only  be  fractal  over  a  certain  range,  so  the 
choice  of  inner  and  outer  scales,  £0  and  £i;  over  which  the  slope  is 
determined  may  be  narrower  than  those  defined  by  purely  objective 
limits . 

Digitized  data  imposes  an  objective  limit  on  el  because  the  curve 
shape  between  discrete  data  points  is  unknown,  so  it  is  simplest  to 
assume  a  straight  line.  As  before,  one  plots  log  N(£)  vs  log£  to 
find  DB.  For  e  <  £m  the  curve  scales  as  a  one  dimensional  line, 
the  assumed  shape  between  adjacent  data  points.  So  the  inner  scale 
of  resolution,  £,,  must  span  at  least  three  data  points,  i.e.,  2£m. 
This  value  happens  to  be  same  as  the  Nykqvist  cutoff  for  Fourier 
spectra,  but  for  different  reasons. 

An  objective  criterion  also  exists  for  the  choice  of  outer  scale, 
£0,  when  the  time  series  contains  waves.  As  shown  in  fig.  2  for  a 
single  sinusoid,  log  L(£)  versus  log  £  shows  nearly  zero  slope  for 
£  <  X/3,  whereas  for  e  slightly  greater  than  X/3  the  plot  is  nearly 
vertical.  So  we  suggest  for  wave  containing  time  series  that  £0  < 
Xm/3,  where  A   is  the  foreseen  minimum  wave  period. 


J.u  q 


4.0  ^ 


"n 1  ',!:':  i  ,!!  I  l 


l  l  :  |  :  i  :  i  I  I  I 

1.5 
I 


n  r 


.0 


2.C 


2.5    3.0 


7  C 

J  .  J 


n  £" 


Figure  2  -  Self-affine  L(e)  for  1800  point  sinusoid. 


Mandelbrot  (1977)  describes  another  way  to  determine  fractal 
dimension  for  a  continuous  curve  in  two  dimensions,  such  as  a 
geographical  coastline.   Here,  curve  length  is  defined  as 


L(e)  =  zN(e) 


(II. 5) 


where  N(£)  is  given  as 


NU)  =  £  ,  (II-6> 


and  F  is  again  the  lacunity. 

Lacunity  measures  the  length  of  one  fractally  scaling  curve  to 
another  of  the  same  D.  I.e.,  assume  that  within  the  limits,  £0  to 
£;  the  fractal  dimension  is  constant  for  any  segment  of  a  real  curve 
of  finite  length.  Then  for  any  two  segments  of  unequal  length, 
II. 6  shows  that  the  step  length  ratio  will  be 


(II. 7) 
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Using  II. 6  in  II. 5,  L(fi)  can  be  written  as 

L(e)  »  Ft1'0  . 


(II. 8) 


Taking  the  log  of  II. 8  yields  the  straight  line  equation, 

log  L(e)  =  logF  +  (l-£>)  log  e  .  (II. 9) 

The  slope  of  log  L(£)  by  log  £  yields  (1  -  D)  .  To  find  D  first 
find  L(£0)  .  Then,  from  II. 9,  log  L(£0)  •=  log  F  +  (1-D)  log£0.  Now 
with  L(£,)  ,  log  L(£j)  =  log  F  +  (l-D)log£i.  This  yields  two 
equations  in  two  unknowns,  D  and  F,  so  in  principle 


log 


LU0) 


D  =   LUl>         .  (11.10) 

log  [-^2] 

This  procedure  amounts  to  spanning  the  curve  with  circles  of 
diameter  £.  So  L(£)  is  also  the  minimum  number  of  circles  of 
diameter  £  required  to  cover  the  curve  times  their  diameter.  Thus, 
L(£)  can  be  measured  by  stepping  along  the  curve  with  a  compass 
whose  span  is  set  to  e.      Mandelbrot  (1985)  terms  this  estimate  the 

8 


"compass  dimension",  Dc. 

To  show  that  Dc  is  equivalent  to  DB,  solve  II. 8  for  D  to  get 

d-  l°*LM  -log'f)  .1 .  Hi. ix, 

log! 

Substituting  L(£)  from  II. 5  gives 

D-    lQ^(£)  -log(F)  (II-12J 

log- 
e 

and  in  the  limit  as  e  -  0,  since  N(e)  is  much  larger  than  F,  this 
yields 

D=   lim  lQg^(£)  ,  (U.13) 

-0   logA 
e 

which  is  DB. 

In  essence  fractal  dimension  is  a  measure  of  the  jaggedness  or 
degree  of  convolution  of  a  curve  given  in  terms  of  the  power  law 
dependence  of  the  curve  length  as  a  function  of  resolution. 

The  box  method  divides  the  data  plane  into  discrete  boxes  of 
length, e,  where  N(£)  is  the  number  of  boxes  having  at  least  one 
point  on  the  curve.  So,  for  DB  every  box  in  the  domain  must  be 
checked,  even  though  most  boxes  are  empty,  while  Dc  considers  only 
the  circles  covering  points  on  the  curve.  So  though  Dc  and  DB  are 
identical,  Dc  is  faster  to  compute. 

C.  Time  series  applications 

Unlike  coastlines,  a  time  series  plot  has  axes  with  different  units 
of  measure,  time  versus  amplitude.  This  poses  a  problem,  since  the 
"length"  of  a  time  series  trace  cannot  be  measured  without 
specifying  an  arbitrary  scaling  ratio  between  axis  units.  I.e., 
different  ratios  will  give  different  values  of  Dc.  Thus,  though 
calculable,  Mandelbrot  (1985)  suggests  that  Dc  may  not  be 
meaningful  for  time  series  analysis,  since  its  value  can  even 
exceed  two,  if  the  y  to  x-axis  scaling  ratio  is  large  enough.  This 
warrants  exploring  a  shift  in  approach  to  measuring  fractal 
dimension,  which  we  describe  as  follows. 

McHardy  and  Czerny  (1987)  applied  a  "self-af f ine" ,  "multi-fractal" 
operator  to  their  time  series  data.   Meaning  "in  the  same  way", 
self-affinity  requires  that  one  axis  be  stretched  or  compressed  at 


each  successive  resolution  by  the  same  constant  factor  with  respect 
to  the  other  axis.  Thus,  McHardy  and  Czerny  defined  their  self- 
affine  "length  metric"  as 


L(t)    =   lfT\F(t  +  c)   -F(t)   |  dt   ,  (11.14) 


e 

o 


where  now  the  abscissa  scaling  changes  by  the  factor,  zn/£n+\,  each 
time  £  changes  by  the  factor  £n+i/£n  (their  length  metric  differs 
from  the  standard  self-affine  length  metric,  which  uses  the 
integrand,  (  [F(t+£)  -F  (t)  ]2  +  £2)  ,/4  ).   This  new  metric  leads  to, 

D   _   dlogL(e)  (11.15) 

d  log  e 


L(£)  from  11.14  will  be  much  longer  at  small  £,  but  note  that  D  in 
11.15  is  defined  as  the  change  rate  of  log  length  with  log 
resolution,  not  the  ratio  of  log  length  to  log  resolution.  Since 
time,  £,  cancels  in  the  length  metric,  L(£)  only  has  units  of 
amplitude.  This  avoids  arbitrary  ratios  between  amplitude  and  time 
units  and  makes  the  problem  one  dimensional,  so  D  is  less  than 
unity.  Thus,  eqn.  11.14  presents  a  natural  choice  for  evaluating 
time  series.   Hereafter,  we  refer  to  D  from  this  method  as   DA. 

L(£)  calculated  from  11.14  will  contain  some  noise  contribution, 
which  is  statistically  independent  of  the  signal.  If  its  level  is 
known,  we  can  estimate  its  contribution  to  L(£)  by  the  formula, 
L2observed  =  Signal  +  ^noise  •  White  noise  has  DA  =  1 ,  so  it  tends  to 
increase  the  fractal  dimension  of  the  time  series. 

D.  Relevant  gravity  wave  and  turbulence  characteristics. 

Gravity  waves  in  the  stable  boundary  layer  can  be  generated  by  a 
number  of  mechanisms,  among  them:  wind  shear  (Kelvin-Helmholtz 
instability),  storm  impulses,  and  flow  over  an  obstacle.  Their 
amplitudes  can  vary  from  a  few  centimeters  to  hundreds  of  meters, 
with  periods  from  less  than  a  minute  up  to  40  minutes  (Stull, 
1988).  Wave  generation  due  to  flow  over  an  obstacle  is  significant 
here.  Hunt  (1980)  gives  the  natural  wavelength  of  flow  over  a  hill 
as 

X    =  2tc!70/BVF  ,  (11.16) 

where  U0  is  mean  wind  speed.  The  hill  "length"  is  given  as  Lu  the 
horizontal  distance  across  the  hill  at  half  its  maximum  elevation. 
If  X  -  5L]  then  lee  waves  may  occur  with  X.  Moreover,  if  e  -  2L,, 
then  strong  lee  waves  are  possible.   Observational  studies  such  as 
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Caughey  and  Readings  (ibid),  Caughey  (ibid),  and  Nai-Ping  (ibid), 
suggest  that  waves  and  turbulence  commonly  coexist  in  the  nocturnal 
boundary  layer,  and  turbulence  generation  can  sometimes  depend  on 
the  wave.  The  idea  that  gravity  waves  transfer  energy  to  Kelvin- 
Helmholtz  waves  which  then  "break"  is  widely  used  as  a  model  for 
turbulence  generation  by  waves  (Stull,  ibid;  Atlas  et  al.,  1970), 
though  this  model  is  not  universally  accepted  (Hines,  1988). 

Finnigan  (ibid)  resolved  the  wave  and  turbulent  portions  of  the 
total  kinetic  energy  by  phase  averaging,  and  showed  horizontal  TKE 
developing  in  the  first  guarter  wave  cycle,  while  vertical  TKE 
arose  during  the  ensuing  third  of  the  cycle.  Both  components 
transferred  energy  from  the  wave  to  turbulence.  They  noted  that 
the  wave  seemed  to  modulate  the  turbulence,  but  the  presence  of  the 
turbulence  had  no  apparent  effect  on  the  wave. 

Gossard  et  al.  (1985)  suggest  a  generation  mechanism  for  boundary 
layer  turbulence  where  local  gradients  of  8,  u,  and  v  increase 
steadily,  while  turbulence  decreases.  Vertical  shear  eventually 
reaches  an  insupportable  value  of  local  Richardson  number,  leading 
to  growth  of  local  Kelvin-Helmholtz  instabilities  which  grow  into 
turbulence.  These  various  proposals  for  turbulence  initiation  and 
wave/turbulence  energy  transfer  are  discussed  in  the  context  of  our 
fractal  analysis  in  section  IV. 
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III.   METHODS 

A.  Data 

The  Boulder  Atmospheric  Observatory  (BAO)  (Kaimal  and  Gaynor,  1983) 
supplied  the  data  for  this  study.  Data  samples  included  eight 
levels  at:  10,  22,  50,  100,  150,  200,  250,  and  300  meters,  of: 

1.   u,  v,  and  w  velocity  components  sampled  at  10  Hz  by 
sonic  anemometers. 

2.  wind  direction  and  speed  at  1  Hz  by  propeller  vane 
anemometers . 

3.  temperature  at  10  Hz  by  platinum  wire  thermometers,  and 
at  1  Hz  by  quartz  thermometers. 

The  1  Hz  data  were  given  only  as  10  second  averages.  The  10  Hz 
data  were  given  at  both  10  Hz  and  as  10  second  averages.  Noise 
levels  were  <  .01  °C  for  the  platinum  wire  thermometer,  and  <  .03 
m  s"1  for  the  sonic  anemometers  (J.  Gaynor,  personal  communication) . 

The  data  and  FORTRAN  input  read  programs  were  sent  on  9  track  tape 
by  John  Gaynor  and  Dave  Welch  of  the  NOAA  Wave  Propagation 
Laboratory.  The  programs  were  run  on  a  NPS  Computer  Science 
Department  Sun4  workstation,  with  analysis  on  the  NPS  mainframe. 
The  data  covered  three  (MST)  time  periods:  A)  2340  9/7/83  -  0340 
9/8/83;  B)   0440  -  0620  9/9/83;  and  C)  0000  -  0620  9/19/86. 

The  temperature  time  series  seemed  the  least  noisy  of  the  available 
data.  So  10  Hz  temperature  data  were  plotted  at  levels  four  and 
five  (100  and  150  meters)  for  all  three  periods.  A  "wave"  was 
apparent  on  visual  scanning  period  A  at  level  four,  between  0040 
and  0105  MST,  September  8,  1983  (see  fig. 3,  windows  21-27). 
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Figure  3  -  Temperature  (°C)  vs.  three  minute  windows. 
The  wave  occurs  in  windows  21  -  27;  turbulence  episodes 
in  windows  45  -  47,  66  -  67,  and  73  -  75. 
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Since  winds  were  southerly  throughout,  the  "wave"  was  consistent 
with  the  above  lee  wave  criteria  and  flow  over  a  small  hill  south 
of  the  tower.  I.e.,  eqn.  11.16  gave  X  —  1900  meters,  for  an 
observed  BVF  =  0.03  Hz,  and  U0  =  9  m  s"1.  From  a  USGS  map  for  the 
120  ft  hill  (elev.  5280  ft  MSL) ,  L{  -  300  meters  along  the  tower 
direction;  so  51^  ~  1500  meters.  This  fits  the  condition  X  ~  5Lj . 
But  before  and  after  the  "wave"  period,  the  BVF  dropped  below  0.01 
Hz,  so  X  exceeded  6200  meters,  well  beyond  the  lee  wave  condition. 
This  account  is  consistent  with  terrain  induction  of  the  wave.  Lee 
wave  conditions  were  approached  for  windows  43-53,  and  61-66,  but 
not  as  closely,  and  waves  were  also  not  evident  in  the  time  series. 
Other  periods  showed  no  clear  wave  evidence.  Period  C  also  showed 
pronounced  noise  spikes  at  20  minute  intervals  and  high  noise 
levels  throughout  most  of  its  length. 

In  period  A  a  reduced  temperature  variance  region  followed  the 
"wave"  episode.  This  let  us  compare  fractal  characteristics  of 
both  "wave"  and  non-wave  portions  of  the  time  series.  Thus,  we 
focused  on  level  four  data  from  period  A  for  this  study. 

B.  Analysis  measures 

Dc,  DA,  turbulent  kinetic  energy  (TKE)  and  its  component  velocity 
variances,  bulk  Richardson  number  (RB) ,  Brunt-Vaisala  frequency 
(BVF) ,  buoyancy  length  (1B)  ,  and  fast  fourier  transform  (FFT) 
spectra  were  all  tested  as  analysis  tools  on  the  u,  v,  w,  and  T 
time  series  at  level  four.  RB,  BVF,  and  1B  are  bulk  measures,  and 
were  taken  across  levels  three  to  five,  a  vertical  span  of  100 
meters.  The  above  parameters  were  checked  for  possible  correlation 
with  DA. 

The  main  disadvantage  of  bulk  measures  is  that  they  apply  to  the 
entire  layer  rather  than  a  local  sensor;  Dc,  DA,  TKE,  FFT  spectra, 
and  variance  are  taken  from  local  values  valid  only  at  the  sensor. 
So  if  an  eddy  is  smaller  than  the  100  meter  layer  thickness,  it  may 
impact  the  local  sensor  but  not  the  heights  at  which  the  bulk 
gradients  of  U,  V,  and  6    are  measured. 

Dc  was  computed  for  the  level  4,  10  Hz  temperatures  for  adjacent 
three  minute  windows  throughout  period  A.  Three  minute  windows 
were  chosen,  since  the  "wave"  had  a  period  of  about  200  seconds. 
As  discussed  above,  et  was  set  to  three  data  points  (.2  sec).  £0 
was  set  to  the  window  length,  1800  data  points,  for  the  first  run. 
For  all  runs  the  time  axis  was  expanded  by  a  factor  of  100,  so  that 
.1  seconds  on  the  ordinate  was  equal  to  .01°C  on  the  abscissa. 

The  first  run  for  Dc  showed  that  the  curve  scaled  linearly  (slope 
on  the  log-log  plots  was  zero)  for  all  log  e  >  2.25,  or  £  >  17.8 
seconds  for  all  windows.  For  most  windows,  the  linearly  scaling 
portion  began  above  e  >  3.2  seconds.  For  better  resolution  on  the 
fractal  part  of  the  curve,  the  run  was  repeated  with  an  s0  =  17.8 
seconds.  Figure  4  shows  a  typical  plot  from  this  run. 
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The  steep,  linear  part  of  the  curve  was  linearly  regressed  to  get 
Dc.  The  linear  region  varied  from  window  to  window,  depending  on 
the  outer  scale.  £0  was  determined  manually  for  each  window, 
according  to  the  following  procedure. 

1.  A  line  was  drawn  with  a  straight  edge  along  the  steep, 
linear  portion  of  the  curve  (point  A  in  Fig.  5) . 

2.  A  second  line  was  drawn  along  the  shallow,  trailing  edge 
of  the  curve  at  high  values  of  log  £. 

3.  The  ordinate  where  these  two  lines  intersected  was  taken 
to  be  the  outer  scale,  e0. 

Some  windows  also  required  that  £,,  be  determined  manually  in  a 
similar  fashion. 
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Figure  4   -  Typical  log-log  plot  from  self-similar 
algorithm  (Dc)  from  the  temperature  time  series. 
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Figure   5   -  Determination  of   largest   £   to  use   in 
calculating  slope  for  Dc. 
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DA  was  also  computed  for  the  temperature  and  the  three  velocity 
components,  using  the  same  three  minute  windows  for  Period  A.  €{ 
and  £0  were  again  .2  and  180  seconds  for  the  first  run.  This 
allowed  approximately  three  decades  of  dynamic  range. 

The  first  run  showed  for  all  windows  that  log  L(£)  vs.  log  e  had 
constant,  non-zero  slope  for  e  from  .2  to  around  60  seconds,  but 
not  beyond  this  time  duration.  Thus,  only  log  e  <  2.75  (s  ~  600) 
was  used,  since  this  still  provided  a  satisfactory  2.5  decades  of 
dynamic  range.  As  mentioned  above,  another  reason  for  setting  £0 
=  600  is  that  Xm/3  =  600  for  the  apparent  "wave"  period.  Thus, 
unlike  Dc,  objective  criteria  could  be  used  to  determine  the  inner 
and  outer  scale  cutoffs  for  the  slopes.  Again  a  simple  linear 
regression  was  used  to  get  DA  for  each  window. 

Mean  turbulent  kinetic  energy  (TKE)  is  defined  as  h{o2  +  a2  +  ow2)  , 
where  au,  av,  and  aw  are  the  respective  standard  deviations  of  u, 
v,  and  w.  Variances  are  often  measured  over  30  minute  to  one  hour 
periods  because  the  apparent  "spectral  gap"  often  occurs  at  about 
1  hour  and  allows  a  convenient  separation  between  large  and  small 
scales.  Statistical  confidence  in  the  value  of  second  moments  such 
as  variance  also  may  not  reach  appropriate  levels  without  averaging 
over  tens  of  minutes.  In  our  case  such  a  long  averaging  time  would 
not  be  useful  in  determining  small  scale  differences  between  waves 
and  true  turbulence.  Instead,  the  TKE  and  variances  were 
calculated  over  the  same  three  minute  intervals  as  Dc  and  DA  for 
purposes  of  direct  comparison.  All  other  parameters  such  as  bulk 
Richardson  number  or  Brunt-Vaisala  frequency  were  also  calculated 
over  three  minute  windows. 

Bulk  Richardson  number,  RB,  was  checked  for  possible  correlations 
with  Dc  or  DA.  Ideally,  Rf/  the  flux  Richardson  number,  or  R,,  the 
gradient  Richardson  number,  should  be  used,  since  our  interest  is 
in  the  local  stability  at  the  sensor.  Both  Rf  and  R;  require  that 
the  local  vertical  gradient  of  the  mean  wind  be  known,  but  this  was 
not  available,  since  the  sensors  were  spaced  50  meters  apart 
vertically.   So,  instead  we  used  the  bulk  Richardson  number, 


Rb   s 


g  \BV  \z 


fl  [(AZf)2  +  (Al/)2] 


where  A  represents  the  quantity  difference  between  the  bottom  and 
top  of  the  layer.  The  layer  thickness,  Az ,  was  100  meters,  with 
the  layer  centered  on  the  level  of  interest. 

The  virtual  potential  temperature,  6V,  was  assumed  equal  to 
potential  temperature,  6.  Though  Stull  (ibid)  emphasizes  that  0V 
can  differ  from  6  by  up  to  4°C,  this  will  not  seriously  affect  RB 
for  three  reasons.  First,  d6v/dz  will  not  differ  much  by 
substituting  9  for  8V;  second,  a  4°C  difference  in  the  denominator 
will  not  change  RB  by  more  than  10  percent;  and  third,  the  absolute 
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value  of 


is  not   important 


In  this   study  RB   is   only 


qualitatively  important  as  an  indicator  of  stability  changes. 

The  Brunt-Vaisala  frequency,  BVF,  was  chosen  to  measure  static 
stability  and,  like  RB,  possible  correlations  with  Dc  or  DA  were 
tested.   BVF  is  defined  as 


BVF  = 


g  ™v 


N  Sv  dz 


As  with  RB,  we  assumed  that  d0v/dz  =  A9/Az  because  we  wished  to 
explore  correlations  in  temporal  changes  in  BVF  with  other 
parameters  and  were  not  interested  in  the  actual  BVF  values. 

1B,  the  buoyancy  length  scale,  was  also  determined  as  the  standard 
deviation  of  the  vertical  velocity  divided  by  the  Brunt-Vaisala 
frequency  ( 1B  =  aw/BVF)  ,  and  was  assumed  to  be  a  measure  of  the 
dominant  eddy  scale. 

Detrended  FFT  spectra  were  also  computed  using  a  cosine  squared,  or 
Hamming  window  to  excise  ringing  and  high  frequency  noise. 
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IV.  RESULTS  AND  DISCUSSION 

A.  Comparisons  with  local  measures 

Figure  3,  the  temperature  series,  shows  both  "turbulence"  and 
"wave"  episodes.  Using  three  minute  windows,  Figure  6a  shows 
distinct  temperature  variance,  aT2,  peaks.  The  largest  occurs  in 
windows  21-28,  cotaneous  with  the  "wave".  Three  other  window 
periods:  45-47,  64-67,  and  73-75  are  cotaneous  with  "turbulence". 
The  other  aT2  peaks,  like  window  62,  were  cotaneous  with  monotonic 
temperature  changes  across  the  window  which  did  not  seem  to  be 
associated  with  the  "wave"  or  "turbulence".  These  are  absent  in 
the  av2  and  aw2  plots  (figs.  6b,  c)  ,  but  do  appear  in  au2  (fig.  6d)  . 
Figure  7  shows  a  TKE  spike  during  the  "wave",  but  also  for  the 
three  "turbulent"  episodes.  So  the  TKE  and  variances  do  not  seem 
to  distinguish  between  "wave"  and  "turbulence";  all  appear  as  local 
maxima. 

For  DA  fig.  8  shows  a  typical  log  L(£)  versus  log  e  plot  for  a 
three  minute  temperature  window.  Most  plots  were  quite  linear  for 
2.5  decades  on  the  ordinate,  from  e  =  .2  to  60  sec.  Figure  9  shows 
the  temperature  DA.  Unlike  variance,  the  lowest  values  occur  in 
windows  21-28  during  the  "wave".  Three  other  local  minima  occur  in 
windows:  42,  63,  and  78,  cotaneous  with  periods  of  seemingly  low 
"turbulence"  in  the  time  series.  Unlike  the  extended  "wave" 
period,  the  other  local  minima  occurred  only  in  isolated  individual 
windows. 

Figure  10  compares  the  temperature  and  velocity  DA.  The  curves 
agree,  with  near-simultaneous  minima  and  maxima  on  all  plots. 
Figure  11  shows  scatter  plots  of  DA  from  one  time  series  versus 
another,  indicating  good  linear  correlation. 
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Figure    6    -   Variances   of:     (a)    T;     (b)    v;     (c)    w;     (d)    u. 
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Figure  7  -  Mean  turbulent  kinetic  energy  (TKE) 
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Figure  8  -  Typical  log-log  plot  from  self-affine 
algorithm  (DA)  from  the  temperature  time  series. 


23 


1.0  -a 


0.9 


0.8 


'0.7  -. 


D 


0.6  -. 


£,0.5  ^ 
D 


0.3  i 


0.2  z. 


0.  1    I  I  I  I  I  I  I  I  I  I  I  T  I  I  I  I  I  I  I  I  |  I  1  I  I  I  I  I  I  I  |  I  I  I  I  I  I  I  |  I  |  II  II  II  I  I  I  |  I  III  III  II  |  II  I  I  I  I  I  I  I  |  I  I  I  I  I  I  I  I  I  | 

0         10        20        30        40        50        60        70        80 

window    number 


Figure  9  -  DA  from  temperature  time  series 
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Figure   10    -   DA   from:     (a)    T;    (b)    w;    (c)    u;    and    (d)    v   time 
series. 
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Figure  11  -  Scattergrams  of  DA  from:  (a)  T  vs.  w;  (b)  u 
vs.  w. 
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Thus,  DA  seems  to  discriminate  between  "turbulence"  (high  DA)  and 
"wave"  (low  DA)  episodes.  For  nocturnal  boundary  layers  this 
points  to  the  possible  use  of  DA  ~  0.35  from  temperature,  or 
DA  «  o.5  from  vertical  velocity  data  as  conditional  sampling  cut- 
off values  to  distinguish  wave  dominated  periods  from  time  series 
of  hot-wire  and  sonic  anemometer  data,  while  retaining  most  of  the 
real  turbulence.  However,  more  data  containing  both  waves  and 
turbulence  must  be  analyzed  before  a  firm  value  can  be  suggested. 
Different  cut-off  values  may  be  appropriate  for  time  series  from 
other  instruments,  such  as  hot  film,  lidar,  and  cup  anemometers 
having  different  amplitude  ranges  and  time  resolutions. 

Though  the  following  suggestion  is  untried  as  yet,  note  that  DA  is 
large  for  "turbulence"  but  small  for  "waves",  while  the  variances 
are  large  for  both.  So  wave  periods  can  be  highlighted  by  scanning 
cotaneous  windows  with  the  combined  operator,  aw:/DA,  while 
turbulence  may  be  highlighted  with  the  operator,  aw2DA. 

All  DA  values  increase  toward  the  end  of  the  "wave",  cotaneous  with 
larger  high  frequency  fluctuations  in  the  T  and  w  time  series. 
This  may  indicate  "wave-break",  or  wave-turbulence  energy  transfer. 
So  DA  may  also  be  useful  in  determining  wave-break  episodes,  unlike 
other  purely  local  measures  such  as  variances  or  TKE.  Waves  may 
augment  the  local  shear,  thereby  instigating  instability,  but  they 
may  also  enhance  TKE  by  increasing  the  buoyancy  flux.  Figure  12 
shows  a  large  increase  in  the  non-phased  averaged  temperature  flux 
and  hence  buoyancy  flux  during  the  latter  part  of  the  "wave" 
period,  without  a  concommitant  increase  in  the  Reynolds  stress 
portion  of  the  TKE  shear  generation  rate.  RB  also  plummets  during 
"wave-break",  and  before  one  "turbulence"  burst,  but  not  prior  to 
the  other  two  "turbulence"  episodes. 

Cotaneous  with  the  calm  preceding  these  "turbulent"  bursts,  two  of 
the  local  DA  minima  occur  just  before  the  "turbulence"  bursts  and 
DA  maxima  at  windows  46  and  66.  This  is  consistent  with  the 
proposal  that  local  shear  may  gradually  increase  during  periods  of 
relative  calm;  but  when  such  shear  reaches  an  insupportable  level, 
K-H  instability  may  be  initiated  and  rapidly  grow  into  turbulence. 
This  shear  based  initiation  differs  from  the  buoyancy  based  "wave" 
case  above. 

The  highest  DA  peak  at  window  36  bears  more  scrutiny,  since  it  is 
cotaneous  with  a  relatively  calm  appearing,  low  variance  portion  of 
the  time  series  following  the  "wave"  period.  One  possibility  is 
that  "wave"  initiated  turbulent  mixing  may  persist  for  some  time  as 
intense  but  small  scale  turbulence,  difficult  to  resolve  by  visual 
scanning  alone.  DA  may  be  a  more  acute  discriminant  than  the  eye 
itself . 
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Fig.  12  -  DaCT),  stress,  and  temperature  flux  in  Period  A. 
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Alternatively,  one  might  also  suggest,  due  to  the  1/e  weighting, 
that  DA  is  largely  a  scale  separator  which  seems  to  distinguish 
between  waves  and  turbulence  only  because  turbulence  in  the  stable 
boundary  layer  is  generally  smaller  scale  than  waves.  However, 
this  does  not  explain  the  DA  minima  also  seen  at  windows:  5,  42, 
56,  63,  and  78,  where  the  fluctuations  seem  to  be  small  scale. 

Since  Dc,  the  self-similar,  fractal  dimension  depends  on  the  y  to 
x-axis  scaling  ratio,  we  arbitrarily  set  .1  second  =  100°C.  As 
expected,  Dc  behaved  differently  from  DA,  rising  dramatically 
during  the  "wave"  episode,  and  falling  immediately  after.  This 
same  behavior  was  observed  during  the  "turbulent"  episodes;  thus, 
Dc  did  not  distinguish  between  the  "wave"  and  turbulence.  Moreover, 
the  long,  ad  hoc  algorithm  involved  a  dynamic  range,  £,  -  e0, 
smaller  and  more  subjective  than  the  limits  used  for  DA.  For  these 
reasons  Dc  does  not  seem  useful  as  a  wave/turbulence  discriminant 
for  time  series  and  was  not  pursued  further. 

Figure  13a, b  show  Fourier  spectra  of  the  temperature  data  before 
and  during  the  wave.  Before  taking  the  FFT,  the  data  were  linearly 
detrended  and  tapered.  Despite  such  processing,  the  plots  are 
quite  noisy,  showing  only  slight  evidence  of  a  spectral  gap  between 
the  "wave"  and  higher  frequency  fluctuations.  There  is  a  decided 
increase  in  slope  for  the  lower  frequency  components  above  0.1  Hz 
which  may  indicate  wave  presence.  However,  the  "turbulence"  periods 
do  not  exhibit  characteristics  of  any  distinctive  sort  in  the  FFT 
plots,  as  shown  in  fig.  14a, b. 

DA  also  has  a  natural  advantage  over  traditional  spectral  methods 
in  that  it  assumes  that  the  amplitude  is  constant  over  the  time 
resolution  in  question.  So  DA  actually  relies  on  local  Chapeau 
basis  functions  rather  than  global  basis  functions,  such  as  Fourier 
series,  Bessel  functions,  etc.  to  obtain  digitized  amplitude 
differences.  Thus,  as  illustrated  in  fig.  15,  DA  has  no  problem  in 
tracking  large,  nearly  discontinuous  jumps  in  amplitude,  while  FFTs 
require  that  linear  combinations  of  wave  numbers  much  higher  than 
£j  be  found  to  both  simulate  such  jumps  and  suppress  fictitious 
ringing  outside  the  area  of  the  jump  discontinuity. 
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Figure  13  -  FFT  spectra  from  temperature  time  series:  (a)  before 
wave;  and  (b)  during  wave. 
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Figure  14  -  FFT  spectra  from  temperature  time  series:  (a) before 
third  turbulent  episode;  and  (b)  during  turbulence. 
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Figure  15  -  Digital  approximation  of  curve  height  at  time 
resolution  e   using  Local  Chapeau  (top  hat)  basis  functions. 
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B.  Comparisons  with  bulk  measures 

Figures  16(a)  and  (b)  show  the  temperature  DA  compared  with  Brunt- 
Vaisala  frequency  (BVF) .  The  BVF  is  highest  during  the  "wave"  when 
DA  is  small.  The  three  DA  minima  at  windows  42,  63,  and  78  occur 
when  the  BVF  is  at  or  near  a  local  maxima,  consistent  with  the  fact 
that  a  higher  BVF  (more  stable)  will  tend  to  dampen  turbulence. 
Also,  most  of  the  DA  maxima  occur  when  the  BVF  is  either  decreasing 
or  at  a  minimum,  which  indicates  a  less  stable,  more  turbulent 
atmosphere.  But  high  BVF  values  occurred  during  both  "wave"  and 
"turbulence"  periods,  so  BVF  cannot  be  used  to  resolve  the  two. 
Figure  16(c)  shows  a  scatter  plot  of  the  temperature  DA  versus  BVF, 
which  indicates  a  weak  negative  correlation. 

Figure  17  shows  another  weak,  but  positive,  correlation  between  DA 
and  1B,  the  buoyancy  length  scale.  The  correlation  is  highest 
during  the  "wave"  episode,  when  1B  drops  rapidly,  simultaneous  with 
a  rapid  drop  in  DA .  After  the  "wave",  DA  and  1B  both  rise.  The  u 
and  T  variances  and  BVF  show  local  maxima  at  window  62,  while  the 
u  and  w  variances  and  1B  show  local  minima.  The  temperature  time 
series  shows  a  rapid  monotonic  temperature  fall.  1B  does  not  seem 
to  distinguish  between  "wave"  and  turbulence.  All  such  episodes 
show  moderate  1B  on  the  order  of  8-12  meters.  Since  the  BVF  is  in 
the  denominator  of  1B,  1B  peaks  during  windows  18-20,  and  35-40  when 
the  BVF  is  low.  This  seems  to  occur  when  the  temperature  time 
series  is  relatively  calm.  Thus,  large  1B  does  not  seem  to 
correlate  with  turbulence,  and  does  not  seem  to  be  a  good  measure 
of  dominant  eddy  scale,  as  some  have  suggested  (Stull,  ibid) . 

Comparison  of  DA  with  RB  in  fig.  18  shows  little  correlation  except 
during  the  "wave"  when  RB  rises  dramatically,  indicating  increased 
stability.  None  of  the  bulk  measures  seem  to  reliably  distinguish 
between  "wave"  and  "turbulence"  episodes.  Note  however,  that  high 
correlations  between  DA  and  BVF,  RB,  or  1B  may  not  be  reasonable, 
since  DA  is  a  local  measure,  whereas  the  latter  three  are  bulk 
measurements  from  widely  vertically  separated  sensors.  Higher 
correlations  might  be  expected,  if  more  local  measurements  of  Ad 
and  Au  were  available.  However,  this  capability  is  not  currently 
present  on  the  BAO  tower  or  similar  facilities. 

The  wave/turbulence  energy  transfer  during  windows  25-28  may 
initially  blend  the  potential  temperatures  across  levels  3  to  5 
(see  Fig.  19),  resulting  in  high  1B  and  low  BVF  after  window  35. 
This  is  consistent  with  the  local  maxima  in  the  w  and  v  variances 
around  window  35.  However,  again  the  large  DA  maxima  in  window  35 
of  the  u,  v,  w,  and  T  traces  may  suggest  continued  turbulent 
blending,  principally  at  smaller  scales,  even  after  the  level  4  and 
5  potential  temperatures  have  become  well-mixed. 

Though  not  suggested  by  our  data,  it  is  still  possible  that  such 
smaller  scale  blending  is  unduly  accentuated,  since  DA  is 
normalized  by  l/£n,  where  n  was  assumed  to  equal  unity.    For 
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inertial  eddy  cascade  processes,  n  may  indeed  be  less  than  unity. 
As  a  first  guess,  since  the  turbulence  energy  spectrum  tends  to 
have  a  -2/3  slope  in  the  inertial  subrange,  one  might  expect  the 
turbulence  velocity  spectrum  to  have  a  -1/3  slope.  In  that  case, 
perhaps  the  correct  self-affinity  factor  for  inertial  subrange 
turbulence  should  be  n  =  1/3.  So  to  modify  DA  one  would  first  have 
to  modify  the  integrand  in  11.14  somehow  to  retain  the  units  of 
pure  amplitude  for  L(e) . 

Secondly,  in  practice  the  inertial  subrange  scaling  assumption  is 
not  entirely  viable.  The  BAO  data  resolution  is  0.1  seconds,  which 
enters  the  dissipation  range  where  the  spectral  energy  slopes  are 
much  higher.  To  eliminate  contributions  from  this  range,  the  inner 
scale  cutoff  would  have  to  be  at  longer  times,  thus  narrowing  the 
dynamic  range  over  which  DA  is  computed.  This  may  be  possible  for 
convective  data,  since  a  longer  window  must  be  used  to  capture 
fluctuations  due  to  boundary  layer  scale  eddies. 

Thirdly,  stability  also  tends  to  compress  the  vertical  extent  of  an 
eddy.  Thus,  inertial  subrange  turbulence  may  not  have  a  -%  slope 
in  the  energy  spectrum,  as  in  self-similar  turbulence.  That  is, 
there  will  be  less  energy  in  the  larger  scales.  Mahrt  and  Gamage 
(1987)  have  studied  vertical/horizontal  velocity  aspect  ratios  for 
structure  functions  as  a  function  of  stability.  However,  their 
stability  categories  were  not  well  defined  by  Rs  or  other  measures. 

In  light  of  these  difficulties  and  lack  of  supporting  data,  we 
retain  the  present  self-affine  operator,  DA,  and  note  that  it 
appears  to  be  the  first  potentially  operationally  useful  wave/ 
turbulence  discriminant.  It  may  also  be  useful  as  a  fast  pre- 
processor to  guickly  locate  intermittent  wave  or  turbulence  periods 
for  more  detailed  phase  averaged  studies. 

As  an  extension  of  this  study,  another  potential  use  for  DA  is  as 
a  general  measure  of  the  degree  of  chaos,  applied  to  systems 
ranging  from  pure  waves  and  limit  cycles,  through  standard  chaotic 
systems  such  as  the  Henon,  Poincare,  and  Lorenz  maps,  as  well  as 
transitional  turbulence,  and  turbulence  displaying  extended 
inertial  subranges.  One  effort  has  begun  to  test  DA  on  a  set  of 
differential  equations  which  describe  coupled  harmonic  oscillators. 
By  eye  the  degree  of  chaos  expressed  in  each  oscillator's  motion 
seems  to  increase  or  decrease  monotonically  along  the  oscillator 
lattice.  An  adequate  measure  this  chaos  has  not  been  shown  as  yet. 
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IV.   SUMMARY  AND  CONCLUSIONS 

Time  series  of  T,  u,  v,  and  w  from  a  stable  atmospheric  boundary 
layer  were  examined  using  fractal  techniques.  The  time  series  were 
taken  from  10Hz  sonic  anemometer  and  hot  wire  data  at  a  height  of 
100  meters  on  the  mast  at  the  Boulder  Atmospheric  Observatory 
(BAO)  .  The  time  series  featured  a  probably  terrain  induced  gravity 
"wave"  period,  as  well  as  several  periods  of  intermittent 
turbulence. 

Both  self-similar  fractal  dimension,  Dc,  and  self-affine  fractal 
dimension,  DA,  were  explored  in  the  analysis.  The  latter  was 
clearly  superior  in  practice  as  well  as  principle.  DA  was  the  only 
measure  among  the  eight  (five  local  and  three  bulk)  studied  which 
seemed  to  distinguish  the  "wave"  from  "turbulent"  episodes, 
reaching  minima  during  "wave"  and  maxima  during  "turbulence". 

Since  DA  is  a  local  measure,  this  implies  that  tower  data  from  a 
single  fast  response  sensor  can  be  used  for  wave/turbulence 
discrimination,  while  phase  averaging  requires  extensive  numerical 
processing  and  microbarograph  array  data  as  well.  Most  tower  sites 
are  not  equipped  with  microbarograph  arrays. 

Phase  averaging  also  requires  that  the  wave  period  be  known  prior 
to  analysis,  otherwise  the  wave  and  turbulent  components  cannot  be 
distinguished.  DA  also  requires  an  estimate  of  the  minimum 
expected  wave  period,  Xm,  to  specify  £0,  the  outer  scale  for  the 
slope.  However,  the  required  precision  is  perhaps  a  factor  of  two, 
much  lower  than  is  needed  for  phase  averaging.  For  Xm  we  note  that 
nocturnal  atmospheric  boundary  layer  gravity  waves  have  typical 
periods  in  the  range  100  -  300  seconds. 

From  the  limited  results  we  tentatively  suggest  a  conditional 
sampling  cutoff  value,  DA  ~  0.35,  to  remove  wave  data  from  a  hot 
wire  temperature  time  series.  However,  more  data  needs  to  be 
analyzed  to  establish  a  firmer  value.  Since  aw2  peaks  during  both 
waves  and  turbulence,  dividing  DA  by  aw2  should  highlight  the  wave 
periods.  Similarly,  turbulence  should  be  highlighted  by  multiplying 
DA  by  aw2. 

Correlations  between  D^  and  various  bulk  measures  of  stability  were 
investigated,  with  weak  correlations  (r  <  .41)  found  with  Brunt- 
Vaisala  frequency  (BVF)  and  buoyancy  length  (1B)  ,  and  little 
correlation  with  bulk  Richardson  number  (RB)  .  The  issue  of 
inherently  low  correlations  between  local  and  bulk  measures  was 
discussed. 

Other  possible  parameters  for  distinguishing  waves  and  turbulence 
were  investigated.  Turbulent  kinetic  energy  (TKE)  proved  less  than 
useful,  since  it  included  wave  energy.  Fast  Fourier  Transform 
spectra  were  able  to  resolve  the  "wave"  event  as  a  change  in  slope 
at  low  frequency,  but  no  distinguishing  features  were  found  for  the 
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"turbulence"  episodes. 

In  terms  of  computational  efficiency  DA  does  not  require  that 
coefficients  be  determined  for  linear  combinations  of  wave 
frequencies  to  approximate  moments  of  the  time  series  amplitude, 
i.e.,  no  transform  is  needed.  Moreover,  since  the  slopes  used  to 
determine  DA  seem  well-behaved,  perhaps  sufficiently  accurate 
values  of  DA  may  be  determined  by  computing  L(e)  only  at  the  inner 
and  outer  time  resolutions.  This  would  render  DA  computations 
orders  of  magnitude  faster  than  FFTs.  Moreover,  FFTs  tend  to  treat 
finite  wave  trains,  i.e.,  discontinuities  quite  poorly.  DA  avoids 
this  problem  because  it  involves  local  Chapeau  basis  functions, 
rather  than  global  transcendental  or  other  basis  functions.  This 
makes  DA  a  potentially  suitable  analysis  tool  for  turbulence 
intermittency  and  coherent  structures  in  general. 

DA  also  shows  promise  in  resolving  wave-turbulence  energy 
transfers,  or  "wave-break"  events.  We  suspect  that  "turbulence" 
following  the  "wave"  period  studied  may  have  been  initiated  by 
growing  positive  buoyancy  flux,  while  other  "turbulence"  episodes 
may  have  relied  more  on  shear  generation  of  TKE.  However, 
correlations  between  DA  and  "wave-break"  should  be  tested  in  more 
detail,  perhaps  with  the  Gossard  (ibid)  data  set  which  involved 
vertical  carriage  traverses  along  the  BAO  tower,  wherein  local 
values  of  Rj  were  measured.  Ultimately,  one  would  like  to  establish 
DA  as  a  function  of  local  stability,  so  that  DA  turbulence  cut-off 
values  could  be  associated  with  critical  Rj. 

Tentatively,  we  interpret  the  anomalous  temperature  DA  peak 
following  the  "wave"  period  as  indicative  of  intense  but  small 
scale  turbulence,  not  easily  resolved  by  visually  scanning  the  time 
series  trace.  Alternatively,  due  to  the  1/e  self-affine 
normalization,  DA  may  act  as  simply  a  "scale"  discriminator  rather 
than  a  wave/turbulence  discriminator.  However,  DA  minima  were  also 
found  during  periods  where  visual  scanning  suggested  low  turbulence 
levels,  not  consistent  with  the  latter  interpretation.  Further 
investigation  to  test  DA  on  convective  data  is  suggested  to  explore 
this  possibility.  Proceeding  in  the  other  direction  as  well,  DA 
should  be  tested  on  systems  which  lie  between  the  wave  and  limit 
cycle  regimes  and  turbulent  systems  displaying  extended  inertial 
subranges,  such  as  the  Henon,  Poincare,  Lorenz,  and  other  standard 
chaotic  attractors,  to  see  if  DA  can  distinguish  degrees  of  chaos 
in  general.  An  immediate  effort  is  underway  to  test  DA  on 
differential  equations  which  describe  coupled  harmonic  oscillators. 
To  the  eye  in  certain  modes  each  oscillator  along  the  lattice  seems 
to  display  increasingly  chaotic  motion.  An  adequate  measure  of  the 
this  chaos  has  not  been  demonstrated  as  yet. 

An  n  =  1/3  self-affine  scaling  factor  for  e  is  discussed  as 
possibly  more  suitable  for  turbulence,  though  this  would  require 
modification  of  the  integrand  in  11.14  to  avoid  making  L(£)  also  a 
function  of  time. 
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